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ABSTRACT 


In this report a frequency-domain Green’s Function Method for unsteady supersonic 
potential flow around complex aircraft configurations is presented. 

We focus here on the supersonic range wherein the linear potential flow assumption is 
valid. In this range the effects of the nonlinear terms in the unsteady supersonic compress- 
ible velocity potential equation are negligible and therefore these terms will be omitted in 
this report. 

The Green’s function method is employed in order to convert the potential-flow dif- 
ferential equation into an integral one. This integral equation is then discretized, through 
standard finite-element technique, to yield a linear algebraic system of equations relat- 
ing the unknown potential to its prescribed co-normalwash (boundary condition) on the 
surface of the aircraft. The arbitrary complex aircraft configuration (e.g., finite-thickness 
wing, wing-body-tail) is discretized into hyperboloidal (twisted quadrilateral) panels. The 
potential and co-normalwash are assumed to vary linearly within each panel. 

The long range goal of our research is to develop a comprehensive theory for unsteady 
supersonic potential aerodynamics which is capable of yielding accurate results even in the 
low supersonic (i.e., high transonic) range. 





LIST OF SYMBOLS 


dm speed of sound 

ai contravariant base vector, see equation (29) 

contravariant base vector, see equation (29) 
0 {M$o - l) 1 / 2 

Djk source integral, equation (40) 

Cjk doublet integral, equation (42) 

E(P) domain function, see equation (19) 

S(p) Dirac delta function 

Sjk Kronecker delta 

finite-element shape function 
G Green’s function 

H Heaviside function 

k reduced frequency, ool/U^ 

£ reference length 

Moo free stream Mach number T/qo / a^o 

n unit normal to 

N unit normal to S g 

N e total number of elements 

N n total number of nodes 

uj circular frequency 

0 non-dimensional frequency, oj£/ 

p point having coordinates x, y, z 

p, control point, (x* , y « , 2 * ) 

P point having coordinates X, Y, Z 

P* control point, (X t ,Y*,Z*) 

p.f. Hadamard finite part 

R P - P t 

R 1 hyperbolic radius, see equation (14) 

t time 

T nondimensional time 

U 0 o velocity of undisturbed flow 

x,y,z space coordinates 

X, Y, Z nondimensional Prandtl-Glauert coordinates 
X = x/0l,Y = y/i,Z = zfl 
a surface of body in x,y,z space 
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£ surface of body in X, Y, Z space 

£, surface of element i in X, Y, Z space 

4> perturbation velocity potential 

<f> nondimensional perturbation velocity 

potential, <f>/U 0 0 t 

$ = (feeder ~Mcx>x) 

■0 co-normalwash in x,y,z space 

$ co-normalwash in X, Y, Z space 



In this report we demonstrate how the Green’s Function Method of Potential Aerody- 
namics may be implemented in the frequency domain so as to enable it to handle imsteady 
supersonic flow around complex aircraft configurations. 

We focus here on the supersonic range wherein the linear potential flow assumption is 
valid. In this range the effects of the nonlinear terms in the unsteady supersonic compress- 
ible velocity potential equation are negligible, and therefore these terms will be omitted in 

this report. 

The Green’s function method (Ref. 1) is employed in order to convert the potential- 
flow differential equation into an integral one. This integral equation is then discretized in 
space, through standard finite-element technique (Refs. 2 and 3), to yield a linear algebraic 
system of equations relating the unknown potential to its prescribed co-normalwash on the 
surface of the aircraft. The arbitrary complex aircraft configuration (e.g., finite-thickness 
wing, wing-body-tail) is discretized into hyperboloidal (twisted quadrilateral) panels. The 
potential and co-normalwash are assumed to vary linearly within each panel. 

1.1 A Brief Description of the Green’s Function Method 

Before getting into the specifics of this report we begin with a brief description of 
the Green’s Function Method. This method applies to the equation of the perturbation 
velocity potential. The potential function $ at any point F* in the flow field is given by an 
integral of terms containing the value of the potential and its co-normal derivative on the 
surface, a , surrounding the body and its wake. An integral equation for the potential on the 
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surface of the body is obtained by letting the point F* t approach a point on the surface. 
With this method, the wake is a natural by-product and is treated as a layer of doublets. 
It may be noted that the integral equation does not require that the boundary condition 
on the co-normalwash be satisfied, but rather makes use of the continuity of the potential 
as the control point approaches the surface o. The tangency boundary conditions are 
automatically satisfied by the type of representation obtained with the Green’s Function 
Method. 

In current applications, the surface of the aircraft is divided into small quadrilateral 
elements. Each element is replaced by a paraboloidal hyperboloid surface defined by the 
four corners of the element. In this process the continuity of the surface is maintained 
but discontinuities in the slopes are introduced. The aircraft wake, on the other hand, is 
divided into strips parallel to the streamlines. These wake strips originate from the trailing 
edge and extend to infinity downstream. It should be noted that integrals over these wake 
strips may be carried out in an analogous way to their subsonic counterpart (see Refs. 5 
and 6). 

In the 0 th order theory, the unknown $ (in the Prandtl-Glauert Space) is assumed 
to be constant within each element, while in the 1st order theory # is taken in the form 
4> = # 0 -f £#1 + r/4 > 2 + £* 7^3 where (£,/?) are local element-wise surface coordinates, and 
the coefficients . . . , $ 3 are chosen to interpolate the $ values at the four corners 
of the element. In either situation the integral equation is approximated by a system 
of algebraic equations. This system of algebraic equations is then solved by standard 
numerical methods. It has been found (see Ref. 4) that in the supersonic range at least a 
1-st order theory is required in order to yield a nonsingular set of algebraic equations due 
to a numerical rather than physical anomaly. 


2. UNSTEADY SUPERSONIC FLOW 

Our point of departure is the linearized equation for the unsteady potential compress 
ibie aerodynamic flow 

1 d?<j) 




ale dt 2 


0 


( 1 ) 


where V„ is the Laplace operator in the physical space while <f) is the perturbation potential. 
Choosing a frame of reference such that the undisturbed flow has velocity 'U in the 
direction of the positive x-axis, the linearized total time derivative is given by 


d (') 


d 


dt~ dt + Uoo dx - 
Introducing the generalized Prandtl-Glauert transformation 

X = x/01, Y = y/l, Z = z/l , T = a^pt/l, $ = cfr/U^l 


(2) 


f The bar is used herein to indicate vector quantities. 


(3) 



where l is a characteristic length, = U^/cioo and 0 = ( M — l) 1 / 2 , Eq. (1) yields* 


VoV$ + ^ ; 


,d 2 $ 


+ 2 M 0 


a 2 $ 


o 


( 4 ) 


ar 2 ' °°dXdT 
where o stands for the supersonic dot-product defined as 

a o b = a x b x - a y b y — a z b z (5) 

where a. and b are two arbitrary vectors. Thus the operator V o V stands for 

a 2 a 2 a 2 


V O V 


( 6 ) 


dX 2 BY 2 dZ 2 

Note that a o a is not necessarily positive. We define the ‘super-norm’ of a vector a as 

<.' = 11*11 = |*o a|'/» (7) 

and will use this notation in later sections of this report. 


2.1 Oscillatory Supersonic Fiow-The Potential Equation 

Since in this report we shall be dealing exclusively with oscillatory flow, it is convenient 
at this point to introduce the complex potential via the equation 

#(x,y,z,r) = $(x,y,z)e* n < r ~ M °°*> (8) 

with 

0 = = kM °°/P ( 0 ) 

U 00 P 

where k — is the reduced frequency, and tn is the circular frequency. 

With this notation, equation (4) may be rewritten asj 

VoV$ + n 2 i = o (io) 

We remark that the case Q = 0 corresponds to steady supersonic flow. 


2.2 Supersonic Integral Equation 

In order to obtain the Supersonic Green’s Function integral equation we proceed as 
follows: With P and P* representing the sending and receiving points respectively, the 
Green’s Function G for Equation (10) satisfies 

Vo VG + n 2 G = 6(?-P„ T-T t ) (11) 

* See Appendix C for derivation. 

t The circumflex is used herein to indicate complex quantities. 

| See Appendix C. 
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G = 0 at oo 

one well known solution of which is given by (see Ref. 4) 

H 


G = 


2 nR' 


cos (UR 1 ) 


where 


Here we define 




if X — X, < 0 and R o R > 0 
otherwise 


it = P - p t = (X - X,)i + {Y- Y t )j +{Z- Z t )k 
R' = PH = l(X - X.) 2 - (Y - Y,)‘‘ - (Z - Z,) 1 1 1/2 


( 12 ) 

(13) 

(14) 


Here {P\H(P, Pi) = 1} defines the zone of influence, or Mach forecone, with vertex at P, . 

Multiplying equation (10) by the Green’s function G and subtracting equation (11) 
multiplied by gives 


g(v o + n 2 6) - $(v ovg + n 2 G) = -s(p - p, 

Making use of the identity 

Vo (aV6) = Vao V6 4- aVoV6 

equation (15) reduces to 

V o (G'V# - *VG) = -6(P - P t )S(P 0 ) 


(15) 


(16) 


(17) 


Next for a closed bounded surface E bounding a volume V, we define the domain 
function 

V 

V 


E(P ) = { 1 if P* $ 
{ * } 10 if P* € 


(18) 


Note that for P* on E the function P(P*) will measure the so-called supersonic solid angle 
of E at P* (see Ref. 4 for details). Hence E(P t ) satisfies the notation 


E{P 4 ) = 1 


where N is the outward unit normal to E. 
It can be shown that 


2; tR>) 


(19) 


JJJf* EdV =&f Nd 2 


( 20 ) 
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for ail / continuous where VE is taken in distribution sense. Multiplying equation (17) 
by the domain function and integrating over the whole space yields 

OO 

JJjEVo (GV# - $VG)dV = -E(P.)&(P t ) (21) 

• ~00 

By suitable integrations by parts equation (21) yields 

OO 

fjj VE o (GV$ - §VG)dV = £(?*)$(?*) (22) 

— OO 

Finally making use of equation (20), equation (22) yields 

E(Pj*(P.) =ff^No(GV* - $VG)dE (23) 

If we utilize equation (12) in (23) we obtain 


2t TE(P t )Q{P t ) = jj ^ cos ClR' No V$dE 

- o V ( cos HP') dS 


(24) 


and the above expression is equivalent to equation (A-12) of Morino (see Ref. 5). 


3. NUMERICAL FORMULATION 

In this section, a space discretization procedure will be introduced in order to approx- 
imate the integral equation by a linear algebraic system of equations. Solving this linear 
algebraic system of equations yields the desired perturbation velocity potential solution on 
the aircraft surface. Once the velocity potential is known, the pressure coefficient may be 
computed through Bernoulli’s Theorem. 


3.1 Finite Element Formulation 


Assuming that the surface E is divided into N e small finite elements E t , Equation (24) 
yields 


2 tt E(P t )${P t ) = Y,Jf |^cos UR' No V$dE 

l ~ 1 


2. 

N, 


u 


$ N O V ( ~ cos UR' ) d£ 
n' 


(25) 


Each surface element E, is approximated by a hyperbolic paraboloid given in the form 

p = Pc + Pi £ + Pin + P-itn (26) 
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where P C ,P] , P 2 and P 3 are obtained in terms of the locations of the four corner points 
as (See Fig. 1) 


Pc \ 


( 1 

1 

1 1 \ 

Pi 

_ 1 

1 

1 

-1 -1 

Pji 

“ 4 

-1 

1 

1 -1 

\Psj 


V-1 

1 

-i i ; 



(27) 


Here, F + _, for instance, refers to the element corner for which £ = +1 and r/ = —1. The 
other corners, P + + , P_ + and P__ are defined similarly. It may be noted that the surface 
defined according to Eq. (26) is continuous since adjacent elements have in common the 
straight line connecting the two common corner points. 


3.2 Surface Geometry for Hyperboloidal Elements 

We note that the geometry of the hyperboloidal element is a particular case of the 
general equation for a surface in three- dimensional Euclidean space which is given by 

P = P(Z,V) (28) 

where £ and r) are generalized curvilinear coordinates on the surface. For a hyperboloidal 
surface, the two basis vectors aj and a 2 are given by 

dp — 

«. = QJ = Pi + ^3 
BP _ __ 

«2 = = P2 + (P 3 (20) 

(See equation (26). 

The unit normal to the surface is given by 

N = <i\ x « 2 /|«i x «-2 1 (30) 

and is directed according to the right hand rule such that the normal points outward from 
the surface (see Fig. 2). The surface element dE is given by 

dE = |djd£ x a 2 d??| 


or 


dE = |ai x a 2 |d£d?? 



3.8 First Order Space Discretization 


In what follows we shall take the potential function $(£, r)) over a surface element, say 


E j, as 

^*(£, h ) = [(1 + £)(1 “ h )$+- + (1 + £)(1 + r /)# + + 

+ (1-0(1 + »?)*-+ + (1 -£)(!-»?)$ . _]/4 


(32) 
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Equation (32) expresses the values of 4>* at any point (£, rj) of E t - in terms of the values of 
# at the four corner points. 

More generally, consider the first-order global shape function with the following defi- 
nition: 


(l + 0(l-//)/4 
(l + £)(! + ri )/ 4 
FikU,rj) =<(1 + £)(l + n )/ 4 
(1 - £)(1 - f?)/4 


if node k coincides with corner 4-— of element i 
if node k coincides with corner ++ of element i 
if node k coincides with corner — + of element i 
if node k coincides with corner — of element i 


otherwise. 


With defined by Equation (33), Equation (32) may be rewritten as 

N n 

**'(£,»?) = E^ti’«)*k (34) 

l 

where N n is the total number of nodes on the surface E, and $ k denotes $ at the k^ 
node. 

Similarly, the supersonic co-normalwash $(= N o V$) may be represented by (See 
Ref. 5) 

&U,r/) = E r/)&t (35) 

k—1 

The same first-order finite-element approximation, Equation (33), has been employed for 
\&. Note that if $ is approximated by the l s ^ order finite-element expression while $ is 
represented by a 0^ order formula, $*(£, rj) — & 2 (0,0) = const., a mixed type formulation 
would result. In subsequent portions of this report, and $ will be taken l s ^ order. 


: Numerical Approximation of the Integral Equation 
Making use of the equations (34) and (35), equation (25) may be rewritten as 

B(p.)i(p.) = EE JL F -t§'™ snn ' dE 

+ E E [" a JL t F »* ° v (§ cos n *') dE . 

N n N n 

E(p.)$(p t ) = e B *^k + E c *k®k 

k — 1 k ~ 1 


1] 


where 

N c 1 .. TT 

= // P^cosQP'dE (37) 

t=i * 

c. fc = - E 2 hfL FikW ° v (w cos ni2 ') ds ( 38 ) 

2=1 4 

are the supersonic source and doublet integrals, respectively. The integrals summed in 
Eqs. (37) and (38) are zero except for those indices i which correspond to elements on 
which node k lies. 

If we now select P* to coincide with the node j so that 

1? = \{X - X 3 ) 2 - {Y - Yj) 2 - {Z - Zj) 2 ^ 2 
then equation (36) yields 




£ + E C )A> 


or in matrix notation 

J {4 t } = [/?,*]{¥*} (39) 

where b^ is the Kronecker delta. 

Equation (39) is the derived numerical approximation of the oscillatory supersonic 
integral equation (24) of Section 2. 

In order to use equation (39) we need to evaluate the following coefficients: 

{ Bjk first order supersonic source coefficients, see equation (37), 

Cjk first order supersonic doublet coefficients, see equation (38) 

By using hyperboidal surface geometry these integrals may be written as 

Bjk = 2 n iff cos nR ' F ik{£> X o, 2 \didri (40) 

2 = 1 * * 


Cjk = “ 5 - E jjf, F ikti>v)N ° v cosfiJS') |gi x a 2 \d£dr) (41) 

Making use of equation (30) the doublet integral (41) becomes 

C J k = -j-Y / JJ^Fik(S,ri)(a j xa 2 )oV(|cosnP')d^ (42) 

a— 1 * 

The appearance of the Heaviside function H under the integral restricts the integration to 
that portion of the panel within the Mach forecone. 
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4. PANEL INTEGRALS FOR THE SOURCE AND DOUBLET COEFFI- 
CIENTS 


The supersonic oscillatory Green’s Function, equation (12), Section 2, involves an 
integrable singularity on the Mach forecone, so that the source integral, equation (40), 
may be interpreted in a classical sense. On the other hand the doublet integral, equation 
(41), involves a derivative of the Green’s Function and one must view it in distribution 
sense in order to properly interpret that integral. 

With this in mind, we shall formulate the procedure for calculating the supersonic 
coefficients. 


4.1 Some Definite Integrals 

In what follows we shall allow the r? integral to be evaluated numerically and will 
analytically compute the £ integral. 

By using equation (29) of Section 3.2, we deduce that 



|«1 X < 22 1 — |(P, + r)P 3 ) x (P 2 + £P3)| 

( 43 ) 

We now 

make the approximation 



M*?) xa 2 (0| ~ Mv) + #(*?)£ 

( 44 ) 

where 


A (ri) = \di{rj) x a 2 (0)| 

(45) 

and 


x d-2 ( 1 ) | ~ M??) X a 2 ( — 1)| 
2 

(46) 


Making use of Equation (46) and recalling the relationship (31), the integrals (40) and 
(41) may be reduced to the consideration of the following ^-integrals: 


V H 

«m(h) = P-f- / C'jftjydi 

-t 1 

0,M = J C s ,di 


(47) 


and 


i i 

iM = / CHR'di 


1 


for m = 0, 1, 2 
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The 0 m {rj) and 7 m (r;) are convergent integrals while the ‘p.f.’ in the ot m (r)) integral 
indicates that these integrals must be interpreted in the sense of the Hadamard Finite Part 
in order to assume a finite value. 


4.2 First Order Source Coefficients 

Note that the function cos QR' / R' is an analytic function of R' except for a pole of 
order 1 at R' = 0. If we represent this function by a Taylor series about zero, we obtain 


• *<*•> 

where h(R') is an analytic function of R' with h(Q) = 0. 

Explicitly 

; _ | ( co » nR ' ~ 1 + V 2 R' 2 /2)/R' for R' ± 0 

l 0 for R' = 0 

Furthermore, it is more convenient to write Fn c (£,ri) as 

Fi k ((, «) = FUM + (Fl(v) 

where 


(48) 


(49) 


(50) 


*a(i) = [^it(i,'i) + ri t (-i,i))]/2 

F, :Ul) = [F, t (l,r,)-F,t(-l, n )}/2 

Then, it is easy to show that first order source integral given by equation (40) may be 

expressed as 

^ ■*! 


B 


1 

> k ~ 2ir 


MvjFikPoWdr) 


t=i 


+ 


-t-i 

+ j ^Ain) + F^B(r,)^(n)dr, 

-l 

) 1 2 +1 
f F ik B {n)fh{v)dri - J A(r))F t Q k ^ 0 (ri)dr) 


n 2 

2 


/ [FlkMri) + F? k B{'n) 7i (v)dr) 

4 1 

f F Ik B {v)l2{n)dr)} + S jk 


n 2 

2 


( 51 ) 
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Thus to evaluate the supersonic source integral we must be able to compute the integrals 
(47) for m = 0,1,2. 

We shall return to the source integral a bit later. 


4.3 Doublet Coefficients 

The integral C,* given in equation (42), Section 3 will generally be singular when we 
deal with panels only partially within the Mach forecone, so that, by a simple calculation 
the Hadarnard finite part of Cjk may be given by 


N e + 1 -+ 1 

p f- 2?E / / (,v) 

*= 1-1 


f'* 


] -1 


cos n R' + UR 1 sin UR' 
R' 3 


R • a i x d 2 (53) 


where R is given by equation (14). 

By using equations (26) and (29) of Section 3 we can write 

R ■ a, x d 2 = M 0 (s?) + CM^ry), (54) 

where 

M 0 {ri) = P 0 • (Pj x P 2 + r\P 3 x P 2 ) 

M l (r,)= P 0 ■ P, x P 3 + rfP 2 • ?i x P 3 
= (Po + r/P 2 ) ■ (P l x P 3 ) 

with 

Po = Pc - P> 

where P c is defined by equation (26). 

We now focus our attention on the function 


/w 


cos Viz -f- Or: sin Viz 


(55) 


This function is analytic except for a pole of order 3 at z — 0. If we represent sin Hz, and 
cos Viz by a Taylor series about zero, then we may write 




n 4 

8 


2 + g{z ) 



where g(z) is an analytic function of a with <?(0) = 0. 



Explicitly 


(cosOz + ftesinriz- 1 - Vi V/2 + flV/8)/z 3 for z ^ 0 

g(z) = { ( 57 ) 

0 for z = 0 


Using the notation 0 m (rj) and ')m(ri) as introduced in equations (47) and also 

introducing 

uM = F? k M 0 (r,) 

Vitln) = FbMoiv) + F&M^n) 
mM = F? k Mi(ri) 

we will be able to write down the doublet integrals Cj k . 

The expression for the doublet coefficient is 


1 T -^ ( f 1 r 1 

= 2 Uik ^ a ° ( r /)^ + J v tk{v)<*i( r l)dTi 

i=l 

f 1 rt 2 r 

+ J x WiM^Wdr) + -y- WiMl3o{rj)dri 

+ f v ik (ri)Pi{ri)dri + [ u> lk (r)){-} 2 {r))dr) 

J \ J — ] 

o 4 +1 ' fl 

- g [/ / Gfc(*?bi(h)eb/ 


+ 1 

+ J WikM^Wdri | + D )k 

-i 


(58) 


The last term on the right hand side of equation (58) is a non singular integral and may 
be integrated numerically. Explicitly we have 


jv e + 1 ■+ 1 

D jk = E / / HF lk tt,rt)R.a } x a 2 g{R')dtdr, 

J ;: 1 - ] -l 


( 59 ) 


From the foregoing it is clear that the task before us is to develop a scheme for evalu- 
ating Bj k and C jk . 

Our approach will be to first compute the « m (??), /? m (i?) and 7 m (r?) for m = 0,1,2 
analytically and then carry out the ^-integrations involving these functions in (51) and 
(58) using a Gaussian quadrature technique. The integral expressions S jk and D jk defined 
by (52) and (59) are to be evaluated numerically in both the f and r/ directions using 
Gaussian quadrature. 
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In the next .section we focus our attention on the indefinite integrals in £ associated 
with o m {ri), 0 m (ri) and <7 m (r?), rn = 0, 1,2. 


5. SOME INDEFINITE INTEGRALS 


In this section we shall explicitly obtain the indefinite integrals 


&*(£,»?) = / r di 

PmttPl) = j £ 


I ?' 3 
IV 


(60) 


and 


'7m(t,r,)= f e n R'd( 


for m = 0, 1, 2 

which occur in the process of evaluating the first order supersonic source and doublet 
coefficients. 

For convenience we write 


R' = \\R\\ = (a^A&C + c ) 1 / 2 


(61) 


where 


a = aj o ai 

b = 2(P 0 + r)P 2 ) o gj 

c = (P 0 + r/P 2 ) 0 (P 0 + r/P 2 ) 


We also let 


d — b 2 — 4ac = — 4(R x dj o jR x Sj ) 
From standard integral tables we obtain 


(62) 


( 63 ) 


o 0 




— (4a£ + 2b)/dR', 
-1/R'{2a£ + 6), 

£ c ~ 3 / 2 


for d ^ 0 

for d = 0 
a > 0 

for d = 0 
a = 0 
c > 0 

otherwise 


(64) 
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— (a/?')" 1 - (^b/2aja 0 

for a 0 

d 9 * 0 


Oil (£,> 7 ) = 

f £dz i 
' J R' 3 " 

2R'/b 2 + ^2c/6 2 )/J?' 
— (gJ ?') -1 — (b/2a)a 0 

for o = 0 
d ^ 0 

for a > 0 
d = 0 

(65) 



£V 3 / 2 /2 

for a = 0 
b = 0 
c > 0 




0 

otherwise 


A third basic integral we shall need 

is 




/n|( 2 /?'-y/® + 2a £ + &)/\A|/\/“ 

for a > 0 
d^O 



— tan ” 1 

2 g£ + b ) /2 yf^aR^J / v / — « 

for a < 0 
d 0 


II 

II 

Sr 

Aaj* 

0 

> 50 , 

/n| 2 ./? , 'y/G + 2 g£ + b\ / 

1 

1 

for g > 0 

d = 0 

(66) 


2R'/b 


for a = 0 
d £ 0 



£/v^ 


for a = 0 
d = Q,c > 0 



0 


for otherwise 



In terms of a 0 , «] and fi a above we are able to express the integrals b' 2 ,/ 6 'i,' 7 0 

for a ^ 0 


«2 (£,*?) = J 


edt 

j ?' 3 


0o - ba i - cot 0 )/a 
(2/3 )i ?' 3 - 4c/*" - 2c 2 / R 1 1 / 6 3 

e 3 c“ 3 /2/ 3 


v 0 


for a — 0 
b 0 

for a = 0 
6 = 0 , c > 0 

otherwise 


(67) 
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!h((,n) - J 


%(£,») = / «'<<£ = 



' 2i?' - ^ 0 6/2« 

for a / 0 


II 

73 ftj 

(2/3)R' z /b 2 — (2/b 2 )R'c 
t 2 c- l l 2 / 2 

for a = 0 
b ^ 0 

for a = 0 
&- 0,c > 0 

(68) 

■ ' ■" I 

■ 0 

otherwise 


' (2a£ + 6)i2'/4a + (4ac - b 2 )0 o /8a 

for a / 0 
d / 0 


(2/3b)(b£ + c) 3 / 2 

' 

0 

II 0 

U'H. 

*2 

(69) 



for a = 0 
6 = 0, c > 0 


■ 0 


otherwise 



Finally in terms of /?„,/?! and 7 0 we obtain 7 j ,/?2 ant * 12 

f [(2/3)#' 3 - *ry 0 ]/2a 
[ (2/5) J2' 5 - ( 2c/3)R ,z }/b 2 

e 2 v^/2 


7i(£,t?) = / 


02 = j 


i l di 
' R! 


{(/2 a)R' - {3b/4a)0 x - (c/2a)/? 0 
(2/5) (6^ + c) 5 / 2 /^ 3 - { 2c / b )0i ~ i c2 / b2 )0o 

e/tsfc 


for a / 0 

for a = 0 
b / 0 

for a = 0 
b - 0,c > 0 

otherwise 

for a /: 0 


(70) 


(71) 


7 2 (C,7) = / 


( (£ - bb/6a)R 13 /4a + (5& 2 - 4ac)i 0 /l6a 2 
[(2/7) J?' 7 - (4/5 )cJ?' 5 + (2 / 3)c 2 R' 3 ] / b 3 

eyfc/3 

0 


for a = 0 
b / 0 

for a =- 0 
6 = 0 
c> 0 

otherwise 

for a / 0 

for a = 0 
b / 0 


for a = 0 
b = Q,c> 0 
otherwise 

( 72 ) 
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A 

In sections 6 and 7 we shall see how the indefinite integrals b m , (3 m and 7 m given 
above are to be utilized in evaluating a. m {rt) ,Pm{v) and 7 m (??) for m = 0, 1 and 2. 

0. EVALUATION OF THE INTEGRALS FOR FULL PANELS 

In the case of a full panel, i.e. one in which {(£5*?)! — 1 < £ < +1,-1 < rj < +1} lies 
entirely within the open Mach forecone = {R\X — X* < 0 and R o R > 0} the evaluation 
of the a m ,/3 m and 7 m follows easily. In this situation these integrals are convergent and 
the Hadamard Finite part is not needed. 

The indefinite integrals a m ,/3 m and 7 m given above are to be utilized in evaluating 
<Xm{v),0m(v) and 7 mW for m = 0,1,2. 

a m{d) — O m (l,?j) b m ( I,!?) 

fimin) = Pm( 1,*?) - /L(“M) ( 73 ) 

lm{n) - 7m(M) “ 

for m — 0, 1, 2 

In this situation the a m {r/), fi m (r)),^ m (ri) are analytic functions of 7 for — 1 < 17 < 4-1 
and the numerical computation of the definite integrals involving these functions and 
appearing in Section 4 may be carried out by a standard numerical integration such as 
Gaussian quadrature. 


7. EVALUATION OF THE INTEGRALS FOR PARTIAL PANELS 

In the situation where a panel lies partially within the Mach forecone the evaluation of 
a m{ r l)y0m( r l) and lm{v) takes a bit more doing. We note that integrals d w , m — 0, 1,2 are 
singular on the Mach cone = {i?|X - A"* <0, R o R — 0}. Thus, for £* such that ]?(£*, 7) 
lies on the Mach cone, we must evaluate b m ( £* , 17 ) in accordance with the Hadamard Finite 
Part, We obtain 

p.f. by,, (£,,/?) = 0, m = 0,1 (74) 

On the contrary, the integrals /3 m and 7 W , rn — 0,1,2, are not singular at (£,, ry) so in 
order to calculate the value of these integrals at that point it is enough to plug that point 
into the expressions for these integrals. There are a few provisos however. A problem 
will occur in calculating /? 0 (£*,f?) from expression (66) in the case a < 0 and d ^ 0, since 
/£'(£*, r l) — 0- However the identity 


PoiLn) = - 


tarT 


/ 2u£ + b \ 
\2yf^dR') 


sin 


2a£ -f b 

\fd 


(75) 


for a <0, d ^ 0 together with the fact that 2a£* + b = ±d and d > 0 for £* with i?(£„, rj) 
on the Mach cone shows that 




IT 


2/-a 


sgn (2a£* + b ) 
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for a, < 0, d ^ 0 


(76) 


In addition, we show, in Appendix A, that 

p.f, &,({„,) = { 5RF* S S" < 2 ^- + 6 > 


for a < 0, d = 0 
otherwise 



A A A 

Since the expressions Pi, $2 an< f 72 are a h given in terms of a 0 ,ai,/3 o and R\ 

there is no difficulty in evaluating these functions at £* with R(^,r/) on the Mach cone 
via use of Equations (68) — (72) where applicable. 


We are now ready to investigate our integrals for a fixed r / with — 1 < rj < 1 . We 
look at the intersection of the interval — 1 < £ < +1 with the Mach forecone {R\R o R > 
0, X — X t < 0}. Four cases may occur: 


(i) The intersection is empty. 

(ii) The intersection is a closed interval [£*, £ u ] with < £ u . 

(See Figure 3.) 

(iii) The intersection is a single point with |£ 0 | = 1 but i?(£ 0 , rj) od\ ^ 0. 
(See Figure 4.) 

(iv) The intersection is a single point £ 0 with 

(a) -1 < £ 0 < 1 and R(£ 0 ,rj) oa^=0 (See Figure 5a.) 

(b) |£ 0 | = 1 and R{( 0 , v) 0 = 0 (See Figure 5b.) 


A point (£„, r/ 0 ) where X - AT <0, R o R = 0 and R o = 0 is called a critical point 
if — 1 < Co < 1, — 1 < 'Ho < 1, i.e., (£ 0 , j? 0 ) is a critical point in Case (iv) above. 

Case (i) 

In this case we define 

a m (r)) = 0, = 0, 7 m (»?) = 0 ( 78 ) 

for m — 0, 1,2 

Case (ii) 

In this case we define 

(x m{ r l) ~ d'rn(Cu) ? ?) ~~ V) 

Pm(v) = ~ PmUhV) ( 79 ) 

7 min) = lmUu,ri) - Imtitin) 
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for m = 0,1,2 



where a m , fi rn and are given as in Section 5. We remark however that if either £ u or 
£/ or both lie on the Mach cone the evaluation of the a rn and /3 m for such £* must follow 
(74), (76) and (77) of this section where applicable. 


Case (iii) 

In this situation, we have a limiting situation where either £/ — ► +1 or f u — > —I. The 
functions d m (£, 77) , /3 m (£, 7?) and 7 m (£, 7/) are continuous at such a point and thus we find 
that 

«m(l) = 0 

Mv) = 0 

and 

7m(»?) = 0 (80) 

for m = 0, 1,2 

Case (iv) - (£ 0 , ho) is a Critical Point 

Assume at first that —1 < rjo < +1. Then the equation R o R — a£ 2 + &£ + c = 0 
possesses a double root at £ 0 = —b/2a with -1 < £„ < 1, provided a < 0. At this point, 
d — (5. (Note: If a = d\ oa, > 0, we cannot have a point on the Mach cone with R 0 gj = 0 
unless |j?| = 0 identically.) Now in accordance with Eq. (A. 19) of Appendix A, we find 
that in this case (a < 0, d = 0) 


while 7 m(Vo) — 0 for m — 0,1,2. 

The ct w (7?), m — 0,1,2 behave in a more complicated manner. These expressions may 
be given in the form 


MVo) = 
P\{Vo) ~ 
(h{Vo) = 


7T 


a 


7T 


a 


7T 


a 


e 


(81) 


where: 


at m {v) = ' g (v) + a S rnS{'1 ~ Vo) for m- 0, 1, 2 


al e Hvo)= 0 

(x[ eg {Vo) = 0 ( 82 ) 


and 


a 


reg 




7T 

{-afT 2 
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and the special distributional contribution to a 0 , eq , and a 2 is given by OmS (ry - ry 0 ) where 

rrCT 


sp 


\R ■ dj x a 2 


(83) 
for m — 0, 1, 2 

In the cases where — 1 < r/ 0 < I and |£ 0 | = 1 the only change with the above is that 
the special contributions c4P, ttj^, and are divided in half. 


Thus 


sp 

Oim 


■n a n 


2 | R • a i x a 2 | 

We have not yet considered the situation where \r) 0 \ = 1. Here if |£ 0 | < 1 we set 

if r) 0 R o s 2 > 0 


(84) 


sp 

Otm 


■£o n n/\ R • X a 2 1 if r/ 0 R. a a 2 < 0 


and if both \ri 0 \ = 1 and |£„| = 1 we set 

0 


sp 

(Xm 


if r) 0 R o a,'i > 0 


lCo* K /\R • aj X a 2 | if f] 0 R o a 2 < 0 


(85) 

for m = 0,1,2 

( 86 ) 


8. RECIPE FOR COMPUTER PROGRAMMING 

In this section we summarize the procedure to be used in implementing the mixed 
analytical-numerical evaluation of the source and doublet coefficients B^ and Cjk as given 
in (40) and (42) respectively. 

Step 1 Check if panel is entirely within the open Mach forecone = {R\X — AC < 0 and 
Ro R > 0). If so compute for each rj, —1 < r) < 1, the o m (r/), (3 rn (r)) and O' m(v) according 
to (73). Go to step S. If the panel is not entirely within the open Mach forecone, go to 
step 2, 

Step 2. For each ry, -1 < ry < +1 classify the intersection of — 1 < £ < -fl with the 
closed Mach forecone = (R\R o R >0, X — X* < 0}. Then compute the a TO (r/), /? m (ry) 
and 7 w (ry) in accordance with items (74)-(86). Go on to Step 3. 

Step 3. Evaluate the ^-integrals involving Oi m (r)), 0 m (r)) and 'y m (r /) for m = 0,1,2 and 
appearing in (40) and (42) by a numerical integration scheme such as Gaussian quadrature. 
In the case where there is a critical point (£ 0 , ry 0 ) within the panel and the a m (rj) have a 
distributional component ot & m interpret the integrals in the form: 

+ 1 

J L(ri)o m {rj)dri 

-j 
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as 

4 I 

J L(r])am g (r))dr) + L(r) 0 )a (r/ 0 ) 

- 1 

where and are as given by items (82)— (86). Go on to step 4. 

Step 4. Finally evaluate the double integrals Sjk and Dj in an entirely numerical way by 
Gaussian quadrature in both £ and r] over the unit interval — 1 < £ < +1, — 1 < ?? < +1. 

Care should be taken on two points: 

(i) The integrand should be set to 0 where R o R < 0. 

(ii) Where R' is sufficiently small then both g{R') and h(R') should be 
replaced by the approximation 


and 


without loss in accuracy. 
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g(R') 

h(R') 


n 6 R ' 3 

144 

n 4 5 6 #' 3 

~ 24 ~ 
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APPENDIX A— CRITICAL POINT AND SPECIAL SINGULARITY 
In this appendix we shall study singular integrals of the form 

p - f . / = F a (r,) (A.l) 

where S(£,rj) is a polynomial in £ of degree < 1. 

We focus on the situation where a = «i o Sj < 0 in the neighborhood of a point 
(£o, ri 0 ), -1 < £ 0 < -f 1, — 1 < r) 0 < +1 where R o ~R = 0, R o a\ = 0 and X - X* < 0. 


We shall find that this integral (A.l) exists only in distribution sense as a function of 
rj and in fact takes the ‘value’ 


F.M 


& ( Co ) Vo) 

R • aj x « 2 1 


i— to 


rc6(ri- ri 0 ) 


(A- 2) 


Proof: We study in detail the prototype situation of Fig. 6 where l oa 2 > 0 at (£ 0 ,?? 0 ) 
and the Mach cone intersects rj = +1 at two £ values both with |£| < 1. The situation 
with R o a 2 < 0 at (£ 0 , rj 0 ) with the cone inverted may be handled by analogy. 


We note that for r/j < rj 0 the cone does not interact the line rj = ifo so that clearly for such 

Vu FsiVi) - 0 . 

Further for > r/ 0 the integral (A.l) has a 3/'2-order singularity at two distinct £ 
values on the line rj = rj X . By definition of the Hadamard Finite Part it follows that 
F a (rj j) = 0 here again. 

In order to demonstrate (A. 2) we therefore need only show that 


r F4v)dri=- 3k^R 

J- 1 \R- a-i x do 


\R- d,x a 2 


In fact if we denote by 


W a (ri) = p.f. / j 


R 13 


d£ dij 


ri ° allowable f 


it is clear from previous remarks about F s (ij) that 

1: W s (ri) will be independent of rj for r/ > rj 0 , and that 


(A. 3) 


(A A) 


for rj > rj 0 


2: 


W. 


+ l 

Vi+i) = f drj 
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so that it, suffices for ris to prove that for any r/ > t} 0 


W 8 ( V ) 


SUcVo)* 


\R-ai x «2|l^ 

Interchanging order of integration in (A. 4) we find that for any rji > q 0 


m 


W s {r)i) = p.f. J j ~Jr 3 — dr} d£ 


allowable { r?I> 

We now denote by V(£,r)i) the integral 

*h 


1 ) = P-f- / dr] 


<lo 


We have previously seen that 


/ V(t,r,)d£ = W a {n) 


W„ 


R oR>0 
X - X, < 0 

is independent of r) for t] > r} 0 . 

Let us now evaluate the integral (A. 7) explicitly. 
We note that for any fixed r) > r) 0 we may write 


« \j{iu ~ £)(£ _ it) 


UU,v) 


(A. 5) 


(A.6) 


(A- 7 ) 

for r/j > Y} 0 


R'- 

in the situation of Fig. 6. Here £((rj) and £ u (r/) are roots of \R\ — 0 with £>(r/) < £„(r/). 
If we now write 

we find that 


W S = W S 


m= L 


«- V((,r,)d( 


k, (,):■’ (A ' 8) 

Next making the change of variable £ = (£ u - £i)p/2 -f (£ u + £e)/2 in (A. 8) we obtain 


W s = W, 


-i l 

M = J 


U\ttu — it)p /_ 2 + _( £u + £t ) / 2, r)]dp 

(1 - p 2)'/2 


(A. 9) 
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From (A. 7) we have seen that for r) > r/ 0 the integral (A. 9) is independent of 17 , Now 
letting 17 l r) 0 we of course have Cu{v) j Co and Ce.{ 17 ) f Co and then obtain that 

1 

W* = U(t o ,r, 0 ) J - _ l p2)lJ - 2 = nU(t 0 ,« 0 ) (A. 10) 

Next we compute U{£ 0 ,r)o) in terms of S{£ o ,i 7 0 ). 

We note that 

|«l| 

so setting £ — Co and using L’Hopital’s Rule as 17 | r ) 0 we obtain after simplification 


rr/t \ r R{CoiVo) 

u ( Co , Vo) = lim — jz-77=g— rr 

n-^no — |ai| (i? oa 2 ) 


and then from (A. 10) we have 


w _ S(£.< y 0 )n 

-\ai\[Roa 2 ) 


(A.ll) 


Now in the situation as pictured in Fig. 6 S oa 2 > 0, In the situation where the Mach 
cone is inverted, R. o a 2 < 0 at ( £ 0 , 17 a ), (see Fig. 7) but r/ 2 < f] 0 in (A.6) causing a double 
sign reversal so that we obtain in either case 


p. f . / 


i ?'~ 3 




5 (Co, *?„) 


u‘i|(i? 0 02) 


t=(o 

r\~r)o 


7T 


Hi - r,„) 


(A.l 2) 


To complete our proof we note from Appendix B that at a point {Co-Vo) with Ro R = 0 
and R 0 d\ = 0 we have that 


i? • dj X d 2 | — j«i j|i? o « 2 


so that finally from (A. 12) 


F,(r/) = p.f. J - 


S(C,r/ ) 


R> 3 




S(£o,»fc>) 


| JR • dj x d 2 


€ = «■> 
f) = f)o 


ttA (/7 - rj 0 ) 


(A. 13) 


(A. 14) 


which completes the proof of (A. 2). 

We are now prepared to relate the foregoing to the evaluation of the o m (? 7 ) for m — 
0, 1,2 at a critical point {Co-, Vo)- F° r simplicity we assume —1 < < 1 and —1 < ?/ 0 < 1 

in our discussion. 

We recall that the tt w (r?) are defined as 


a m {ri) = 



for rn = 0, 1, 2 
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In the cases m — 0 or m = 1, £ w will of course be a polynomial of degree < 1 so that 


from item (A. 2) we can immediately say that with S(^ 0 ,rj 0 ) = £ 


m 

o 


ami'll) 


\R-aiX 


Co 

-r) 0 


tin - Vo) 


In the notation of Section 7 we set «m®(/?o) — 0 and 

~nCo n 


sp 

Otrn 


I R - “i x a *\ V't 


for m = 0 or I 


for m = 0 or 1 


(A. 15) 


In the case rn = 2 we are dealing with 


<*M = P-f. I |,3^e 


This integral contains both a regular and a singular component. To isolate them let us 
note that if R' = (a£ 2 + b£ ■+ c) 1 / 2 then £ 2 = R !2 /a + S(£, rj) where S(£,ri) is a polynomial 
in £ of degree < 1. Note that £ 2 = S{£ 0 , r/ 0 ). 


Thus we may write 

= a / W d( + P ' f - / 

Now the second integral on the right of (A. 16) will equal 

-Stto>if 0 )ir I 


(A. 16) 


i? • O] x a'2 


£ ~6o 


% - h 0 ) 


by (A. 2) but since £; = S(£,,ho) this integral must equal 


-ire 

\R ■ d\ x a,2 1 


(-(o 

ri-rjo 


{A. 17 ) 


and we set 

as we have stated in Eq. (83). 
To proceed we must have 


of = . 


jrji 

R ■ a i x d 2 j 


n = r)o 


= 1 / > • 
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But this integral equals 


a 


PoM 


in our notation of Section 7. 
To show 

it stiffic.es to show that 


4 eg (vo) = 


Po{Vo) = 


7T 


(-a) 3 / 2 

7T 

\T^a 


We proceed with the latter. In fact we compute firniw) f° r m = 0, 1,2 at once. Suppose 
the Mach cone cuts our panel as in Fig. 6 and we consider p > r/ 0 . Then with the notation 
R' = yf{tZ u -£)(£- Cl) we see that 


Prn(v) = - 7 = f 

V~ a l 
Ki 


e 


\f{tu - C){C - Ci) 


dC 


(A. 18) 


for rn = 0,1,2 

The change of variable C = (Cu + Cl)/ 2 + (£ u - Ct)p/ 2 transforms A. 18 into 

+i 


Pm{n) = 7 j== J 


i(e„ + f ( )/2 + Uu - (i)p/2] m 


dp for m = 0, 1, 2 


\jl- P 2 

Now letting r/ [ r/ 0 we have that Cu 1 Co and Ct 1 Co so that in the limit 


pm V 
Pm(r)o) = -7= / ; 


dp 7r£ 0 m 


i 0T7 2 ^ 

An almost identical argument shows that 

IrnW = j C m \f(Cu~C)(C~Cl)dC 

(e 


{A. 19) 
for m = 0,1,2 


for m = 0, 1, 2 


+ l 


transforms into 

lrn{v) = \f~(i f [(£„ + Ct )/ 2 + (£„ - 6)p/2p [(e„ - £/)/2] 2 /l - P 2 dp 
■1 

for m = 0, 1, 2 


so that as r/ j r/ 0 , J 7,, 0 T £o we obtain 

lm(Vo) - 0 form =0,1,2 


(A. 20) 
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In this brief appendix we prove an elementary lemma concerning the superdot product. 
This lemma comes into play in proving formula (A. 2) of Appendix A. 

LEMMA: Let a, b and c be three vectors in R 3 . Assume that (i) aoa = Q. (ii) a o b = 0 
and (iii) bob < Q. Then 

|a • b x c| = \[—b o b | a o g|. (U.l) 

We refer the reader to Ref. 3 for definition and properties of the superdot product. 

PROOF : Without loss of generality we may assume coordinates have been rotated so that 
a z — 0. We may assume |a| 0. 

Then a takes the form 

a = a x i + a tJ j with a\ — a'y ^ 0 from (i). 

Let us write 

l = b x i + b y j + b z k 

then 0 = a x b x - a y b y from (ii) or equivalently 

a x h x = a y b y (B. 2) 

We now proceed by cases. 

Case 1: If a x = a y / 0 then from (B.2) it follows that b x = b y . Then 

/ d X (ly 0 \ 

a • b X c = Det J b x b y b z j 

V c x c y c z J 

Q-xiPy^Z b Z Cy) — Q-y(b X C z b z C x ) 

— {^x^x Q-yCy)bz 

Thus a • b X c = (a o c)b z and |a • b x c| = | a o c\\[—b o b in this case. 

Case 2: Here a x = ~a y ^ 0. 

It follows from (B.2) that b x = — b y 
Then again 

/ a x (l y 0 \ 

a • 6 X c = Det b x b y b z I 

\C X Cy C g J 


30 



{by^Z b Z Cy} Gy(6 2 C 2 b z C x ) 

( G z c z — a y c y) ~~ ^2 

(a o c)(—b z ) 


and so 


la • b X cl = |'a o c 


\\/—b o 


b once again. 


APPENDIX C DERIVATION OF THE SUPERSONIC 

OSCILLATORY P.D.E. 

We begin with the linearized potential flow equation: 


V 2 </> 


1 (B rr B \ 2 J /4 

, fu + U oo a_ ) <t> - () 
oo 


aLKBt ' '^BxJ 


(C.l) 


Passing to Prandtl-Glauert coordinates after the introduction of scaled variables as 
indicated in equation (3) we proceed as follows: 


2 1 2U QO B 2 4> Uld*<f>_ 

09 alo Bt 2 alo BtBx at Bx 2 




3 2 <f> -I- 1 <9 2 <^> 2Moo B 2 (j) 


Bx 2 By 2 Bz 2 at Bt 2 a ^ BtBx 


P 2 -r 


B 2 $ 


s 2 $ a 2 # 

+ TiTT-irvi + 


B{PIX) 2 ' 3(IY) 2 B(IZ ) 2 


<9 2 <I> 


2M, 


OO 


B 2 ® 


a 2 o B(lT/ aoo fi) 2 a TO B(IT /a^f^B^lX) 


which results after simplification in 


<9 2 $ 5 2 $ d 2 $ 2 <9 2 $ o d 2 $ 

BX 2 + BY 2 + BZ 2 “ • dT* " 2Moo bxbt 


(C.l) 


At this point we introduce $(X, Y, Z) via the equation 


$(X,Y,Z,T) = S(X,Y,Z)e^ r M °° x > 


3.1 



The consequences of this transformation are as follows: 


OIL - (S® _ iUM $\ e <n(r-M»X) 

ax \ax 00 ) 

S = (lt& ” 2 *nM 00 |J - 

d 2 $ _ d 2 $ iU(T-MocX) 

QY'l QY 2 

® ' 2 ® - d 2 ^ iCl(T- Moo X) 

dZ 2 dZ 2 

df 
d 2 § 
df 2 


in&eWT-MooX) 

~ Moo X) 


d2 l ~ - (n 2 M $ + ,*n— V np, -- M <»*> 
<yxd:r V M °o® + * u dxr 


Then equat ion (C.2) can be written 


d 2 $ 

ax 2 




+ 2f0M 0O ^-y + U 2 Mlo^ + 


d 2 4 d 2 d 

+ 77 


+/? 2 f2 2 4 - 2 H 2 M£>4 - 2iUM 00 j^ = 0 


These equations finally lead to: 


VoV& + n 2 4 = o 


(C.3) 


which is Eq. (10) of the text. 
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